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ABSTRACT 

Wc introduce a class of eccentric discs with "strong" density cusps whose potentials 
are of Stackel form in elliptic coordinates. Our models exhibit some striking features: 
sufficiently close to the location of the cusp, the potential and surface density distri- 
bution diverge as <& oc r^^ and E oc r^^, respectively. As we move outward from the 
centre, the model takes a non-axisymmetric, lopsided structure. In the limit, when r 
tends to infinity, the isocontours of $ and E become spherically symmetric. It is shown 
that the configuration space is occupied by three families of regular orbits: eccentric 
butterfly, aligned loop and horseshoe orbits. These orbits are properly aligned with the 
surface density distribution and can be used to construct self-consistent equilibrium 
states. 
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1 INTRODUCTION 

High resolution observations based on Hubble Space Tele- 
scope photometry of nearby galaxies have increased our un- 
derstanding of the central regions of elliptical and spiral 
galaxies. It was found that in most galaxies density di- 
verges toward the centre in a power-law cusp. In the pres- 
ence of a cusp, regular box orbits are destroyed and re- 
placed by chaotic orbits (Gerhard & Binney 1985). Through 
a fast mixing phenomenon, stochastic orbits cause the or- 
bital structure to become axisymmetric at least near the 
centre (Merritt & Valluri 1996). These results are confirmed 
by the findings of Zhao et al. (1999, hereafter Z99). Their 
study reveals that highly non-axisymmetric, scale-free mass 
models can not be constructed self-consistently. Among the 
models studied for self-consistency, one can refer to the in- 
tegrable, cuspy models of Sridhar & Touma (1997, hereafter 
ST97). Without a nuclear black hole (BH), centrophobic 
bananas are the only family of orbits presenting in ST97 
discs. Although such orbits elongate in the same direction 
as density profile, the orbital angular momentum takes a lo- 
cal minimum somewhere rather than the major axis where 
the surface density has a maximum. This is the main obsta- 
cle for building self-consistent equilibria by regular bananas 
(Syer & Zhao 1998; Z99). A similar situation occurs for anti- 
aligned tube and high resonance orbits for which one could 
not be able to fit the curvatures of orbits and surface den- 
sity distribution near the major axis (Z99). According to 
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the results of Miralda-Escude & Schwarzschild (1989), it is 
only possible to construct self-consistent models by certain 
families of fish orbits. 

The orbital structure of stellar systems is enhanced by 
central BHs in a difi^erent manner. Although nuclear BHs 
destroy box orbits, they enforce some degree of regularity 
in both centred and eccentric discs (Sridhar & Touma 1999, 
hereafter ST99). In systems with analytical cores and central 
BHs, a family of long-axis tube orbits can help the host 
galaxy to maintain its non-axisymmetric structure within 
the BH sphere of influence (Jalali 1999). 

In this paper, we present a class of non-scale- free, lop- 
sided discs, which display a collection of properties ex- 
pected in self-consistent non-axisymmetric cuspy systems. 
Our models are of Stackel form in elliptic coordinates (e.g., 
Binney & Tremaine 1987) for which the Hamilton- Jacobi 
equation separates and stellar orbits are regular. In central 
regions where the effect of the cusp dominates, the potential 
functions of our distributed mass models are proportional 
to r~^ as r ^ 0. So, we attain an axisymmetric structure 
near the centre which is consistent with the predicted na- 
ture of density cusps. The slope of potential function changes 
sign as we depart from the centre and our model galaxies 
considerably become non-axisymmetric. Non-axisymmetric 
structure is supported by a family of eccentric loop orbits, 
which are aligned with the lopsidedness. Our potential func- 
tions have a local minimum around of which a family of 
eccentric butterfly orbits emerges. Close to the centre, loop 
orbits break down and give birth to a new family of orbits, 
horseshoe orbits. Stars moving in horseshoes lose their ki- 
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netic energy as they approach to the centre and contribute 
a large amount of mass to form a cusp. Our models can be 
applied to the study of dynamics in systems with double 
nucleus such as M31 (Tremaine 1995, hereafter T95) and 
NGC4486B (Lauer et al. 1996). 



2 THE MODEL 

Consider the Hamiltonian function 



(1) 



which is described in cartesian coordinates, {x,y). The vari- 
ables px and Py denote the momenta conjugate to a; and y, 
respectively. <E> is the potential due to the self-gravity of the 
disc. Let us express Ti. in elliptic coordinates, (u, v), through 
the following transformations 



X = a(l + coshii cos ii) 
y = asinhusin«, 
u > 0, < < 27r, 



(2) 
(3) 



where a is constant and 2a is the distance between the foci 
of confocal ellipses and hyperbolas defined by the curves of 
constant u and v, respectively. In the new coordinates, the 
Hamiltonian function becomes 
1 



n = 



2a-^(sinh u + sin v 



■{pI +pI) +Hu,v), 



(4) 



with pu and pv being the new canonical momenta. We think 
of those potentials which take Stackel form in elliptic coor- 
dinates. The most general potential of Stackel form is 



F{u) + G(v) 
2a? (sinh^ u ~\- siv? v) 



(5) 



where F and G are arbitrary functions of their arguments. 
By this assumption, the Hamilton- Jacobi equation separates 
and results in the second integral of motion, 72. We get 

l2=pl~2a^Esmh'^u + F{u), (6) 

or equivalently 

~l2^pl-2a^Esmv + G{v), (7) 

where E is the total energy of the system, E = 71. 
We now introduce a class of potentials with 

F{u) = C(coshu)^, 

G{v) = ~Ccosv\cosv['-^, (8) 

where C > and 7 are constant parameters. One can readily 
verify that 

coshii = — (r + s), cosv — —(r — s), (9) 
where 

= + y^, — {x — 2a)^ -I- y^ . (10) 

We substitute from ( |lo|) into (^) and express $ in the [x, y) 
coordinates: 



-I- 

K 



K 



{r + sy - {r - s)\r - s\''~^ 



2rs 



(11) 



C(2a) 



The surface density distribution, associated with <1>, is de- 
termined as (see Binney & Tremaine 1987): 

{V^^)dxdy 



E{x,y') = 



47r2G 



^{x' -xy + {y' -yy 



(12) 



We examine the characteristics of the potential and surface 
density functions for small and large radii. Very close to the 
centre, we have r <C s that simplifies (^l|) as follows 

Ks^-' (l + i)^ + (l- i)^ 



One can expand (1 -I- -)'' and (1 
obtain 



(13) 

in terms of r/s to 



$ : 



Ks'<- 



1 + 



r(7 + i) 



^ (2n)!r(7- 

n— 1 



2n + l) 



(14) 



where F is the well known Gamma function. As r tends 
to zero, s is approximated by 2a and r/s — > 0. Therefore, 
Equation ( [l4| ) reads 

K{2a)''-^ 



<!> 



(15) 



Dimensional considerations show that the surface density E 
will approximately be proportional to r~^. Thus, sufficiently 
close to the centre, we obtain a strong density cusp with 
spherical symmetry. When r tends to infinity, the potential 
$ is approximated as 

/-2 



$ « K2^- 
So, we find out 



(16) 



(17) 



We have to select those values of 7 for which the surface 
density distribution is plausible and orbits are bounded. Ac- 
cording to (0), the surface density decays outward (r 00) 
for 7 < 3. Moreover, Equation ( |l6| ) shows that orbits will be 
escaping if 7 < 2. To verify this, consider the force exerted 
on a star, which is equal to — V$. This force will always be 
directed outward for 7 < 2 and results in escaping motions. 
Therefore, we are confined to 2 < 7 < 3. 

We have used Equations (^) and (^|) to compute "I> 
(Figure 1) and E (Figure 2) for 7 — 2.8. Due to the com- 
plexity of V^$, we have utilized a numerical scheme to eval- 
uate the double integral of (^2|). The potential and surface 
density functions are symmetric with respect to the x-axis 
and are cuspy at {x — 0,y = 0). The potential "1> has a local 
minimum at {x = a,y — 0) that plays an important role in 
the evolution of orbits. This minimum point has no image 
in the plane of the surface density isocontours. The surface 
density monotonically decreases outward from the centre. 
As it is evident from Figure 2, a non-axisymmetric, lopsided 
structure is present at moderate distances from the centre. 



3 ORBITS 

To this end, we classify orbit families. Having the two iso- 
lating integrals E and I2, one can find the possible regions 
of motion by employing the positiveness of and pi in (^ 
and (^). We define the following functions: 



f{u) = -2a'^ E sinh^ u + F(u), 
g{v) = -2a^Esm'^ V + G{v), 



(18) 
(19) 
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(a) 




Figure 1. As a demonstrating example, we liave determined the 
potential $ for 7 = 2.8, a = 0.5 and K = 0.2. Figures (a) and (b) 
show the isocontours and a three dimensional view of <&, respec- 
tively. 



where F{u) and G{v) are given as (H). By virtue of > 



and p„ > one can write 



-h 



fin) 



> 
> 



(20) 
(21) 



Due to the nature of <I>, no motion exists for negative ener- 
gies. Hence, E can only take positive values, E > 0. Our 
classification is based on the behavior of f{u) and g{v). 
The most general form of f{u) is attained for 7C < ia^E. 
In such a circumstance, f{u) takes a local maximum at 
u — 0, fM ~ /(O) = C, and a global minimum at it = itm, 
/m = f{y-m), where 



?^ - 




Figure 2. The surface density isocontours corresponding to the 
potential function of Figure 1. 



coshUiT, — 



and 



C7 



/ni = — 2a^iJsinh^ itm + C(coshMm)^ 



According to (^ 

h > fm. 



we obtain 



(22) 



(23) 



(24) 



On the other hand, g{v) has a global maximum at t; = tt, 
gm ~ gin) — C, and two global minima at w = 7r/2 
and V = 3n/2, gm=g{n/2)=g{3n/2)=-2a^E. Therefore, In- 
equality (Ell) implies 



I2 < 2a E. 

By combining ( [24[ ) and (^ 
/m < /2 < 2a^E. 



(25) 



one achieves 



(26) 

It should be noted that 2a^E > C. This is because of 2 < 
7 < 3. /m and in consequence I2, can take both positive and 
negative values. Depending on the value of I2, three general 
types of orbits are generated: 

(i) Eccentric Butterflies. For C < I2 < 2a E, the al- 
lowed values for u and v are 



1 <V < Vt^2, Vb,3 <V< Vb^i, 



(27) 



where uo and vt^i {i = 1,2,3,4) are the roots of f{u) = 
I2 and g{v) = —I2, respectively. As Figure 3a shows, the 
horizontal line that indicates the level of I2, intersects the 
graph of f{u) at one point, which specifies the value of uq. 
The line corresponding to the level of ~l2 intersects g{v) 
at four points that give the values of v^^iS (Figure 3b). In 
this case the motion takes place in a region bounded by 
the coordinate curves it = no and v = Vb,i. The orbits fill 
the shaded region of Figure 4a. These are butterfly orbits 
(de Zeeuw 1985) displaced from the centre. We call them 
eccentric butterfly orbits. 
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(e) if) 

Figure 3. The graphs of f(u) and g{v) for 7 = 2.8, C = 0.2, 
i? = 1.2 and a = 0.5. The horizontal lines indicate the levels of I2 
and —I2 in the graphs of f{u) and g{v), respectively, (a) I2 = 0.4 
(b) ~h = -0.4 (c) h = -0.25 (d) -h = 0.25 (e) h = -0.1 (f) 
-I2 = 0.1. 



(ii) Aligned Loops. We now let I2 be negative so that 
/m < I2 < —C. In this case the equation /(it) = I2 has two 
roots, iii,i and it;, 2, which can be identified by the intersec- 
tions of /(it) and the level line of I2 (see Figure 3c). The 
equation g{v) = —I2 has no real roots and Inequality (|2l]) 
is always satisfied (Figure 3d). The allowed ranges of it and 
II will be 

ui.i <u< UL2, 0<v <2tt. (28) 




(c) (d) 

Figure 4. The possible families of orbits: (a) an eccentric but- 
terfly orbit (b) an aligned loop orbit (c) a horseshoe orbit (d) a 
lens orbit associated with I2 = C. 

The orbits fill a tubular region as shown in Figure 4b. These 
orbits are bound to the curves of it = it;,i and it — it;, 2 and 
elongate in the same direction as lopsidedness. Following 
ST99, they are called aligned loops. 

(iii) Horseshoes. For — C < /2 < C, we have a different 
story. In this case, both of the equations /(it) — I2 and 
g{v) = —I2 have two roots. We denote these roots by it = 
Uh,i and i; = 11^,1 (i = 1, 2). In other words, the level lines of 
zb/2 intersect the graphs of /(it) and g{v) at two points as 
shown in Figures 3e and 3f. The orbits fill the shaded region 
of Figure 4c, which looks like a horseshoe. We call these 
horseshoe orbits. The orbital angular momentum of stars 
moving in horseshoes (G — xpy — ypx) flips sign when stars 
arrive at one of the coordinate curves v = Vh,i or v — Vh,2- 

For 7C > 4a^E, /(it) is a monotonically increasing 
function of u and eccentric butterflies are the only existing 
family of orbits. There are three transitional cases corre- 
sponding to I2 = C, I2 ~ 2a? E and I2 ~ fm. For I2 = C, 
eccentric butterflies extend to a lens orbit as shown in Fig- 
ure 4d. For I2 = 2a? E, stars undergo a rectilinear motion 
on the line x = a with the amplitude of zbasinhito in the 
y-direction. For I2 = /m, loop orbits are squeezed to an 
elliptical orbit defined by it = itm- 



4 DISCUSSIONS 

In this work we explore a credible model based on the self- 
gravity of stellar discs to explain how an eccentric disc, with 
strong density cusp, can be in equilibrium. Our mass models 
exhibit most features of eccentric stellar systems, especially, 
double nucleus ones such as M31 and NGC4486B. 

All of the orbits of our model discs are non-chaotic. 
Below, we clarify how the existing families of orbits help the 
eccentric disc to maintain the assumed structure. 

The force exerted on a star is equal to — V$. The mo- 
tion under the influence of this force can be tracked on the 
potential hill of Figure lb. This helps us to better imagine 
the motion trajectories. 
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As Figure lb shows, the potential lunction is concave. 
A test particle released from distant regions with x > and 
"small" initial velocity, slides down on the potential hill and 
moves toward the local minimum a,t {x — a,y — 0). After 
passing through the neighborhood of this point (there are 
some trajectories that exactly visit the minimum point), the 
test particle climbs on the potential hill until its potential 
energy becomes maximum. Then, the particle begins to slip 
down again. This process is repeated and the trajectory of 
the particle fills an eccentric butterfly orbit. Stars moving 
in eccentric butterflies form a local group in the vicinity of 
{x = a,y = 0). The accumulation of stars around this local 
minimum of $ can create a second nucleus like P2 in M31 
(see T95). The predicted second nucleus will approximately 
be located at the "centre" of loop orbits while the brighter 
nucleus (PI) is at the location of the cusp. 

Aligned loop orbits occur when the orbital angular mo- 
mentum is high enough to prevent the test particle to slide 
down on the potential hill. The boundaries of loop orbits are 
defined by the ellipses u = uli and u = ui^2- The central 
cusp is located at one of the foci of these ellipses. Aligned 
loops have the same orientation as the surface density iso- 
contours (compare Figures 2 and 4b). Thus, according to 
the results of Z99, it is possible to construct a self-consistent 
model using aligned loop orbits. 

Similarly, we can describe the behavior of horseshoe or- 
bits. Stars that start their motion sufficiently close to the 
centre, are repelled from the centre because the force vector 
is not directed inward in this region. As they move outward, 
their orbits are bent and cross the x-axis with non-zero an- 
gular momentum. These stars considerably lose their kinetic 
energy as they approach the centre (this is equivalent to 
their climb on the cuspy region of the potential hill). Mean- 
while, the orbital angular momentum takes a minimum and 
switches sign somewhere on the boundary of horseshoe orbit. 
This boundary is defined hy v — v^.i (or v — Vh.2) and can 
be chosen arbitrarily close to the centre. These stars spend 
much time near the centre and deposit a large amount of 
mass, which generates a cusp. Therefore, horseshoe orbits 
can be used to construct a self-consistent strong cusp. The 
method of Z99 is no longer applicable to horseshoes because 
such orbits don't cross the long axis (here the x-axis) near 
the centre. In fact, horseshoe orbits are an especial class of 
boxlets that appropriately bend toward the centre. The lack 
of such a property in banana orbits causes the ST97 discs 
to be non-self-consistent. 

In the case of M31 and NGC4486B, if we suppose that 
loop and high-energy butterfly orbits control the overall 
shape of outer regions, horseshoe orbits together with low- 
energy butterfiies (small- amplitude liberations around the 
local minimum of $) can support the existence and stabil- 
ity of a double nucleus. The parameter a will indicate the 
distance between PI and P2. 

There remains an important question: what does hap- 
pen to a star just at the centre? The centre of the model, 
where the cusp has been located, is inherently unstable. 
With a small disturbance, stars located at {x = 0,y = 0) 
are repelled from the centre. But, the time that they spend 
near the centre will be much longer than that of distant re- 
gions when they move in horseshoes. We remark that the 
stars of central regions live in horseshoe orbits. Although 
one can place a point mass (black hole) at the centre with- 



out altering the Stackel nature of the potential, such a point 
mass will not remain in equilibrium and leaves the centre. 
Based on the results of this paper, we conjecture that there 
may not be any mass concentration just at the centre of 
cuspy galaxies. However, a very dense region exists arbitrar- 
ily close to the centre! This may be an explanation of dark 
objects at the centre of cuspy galaxies. The centre of our 
model galaxies is unreachable. Our next goal is to apply the 
method of Schwarzschild (1979,1993) for the investigation of 
self-consistency. 
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